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Abstract. Given an autonomous system of ordinary differential equations (ODE), we consider 
developing practical models for the deterministic, slow/coarse behavior of the ODE system. 
Two types of coarse variables are considered. The first type consists of running finite time av- 
erages of phase functions. Approaches to construct the coarse evolution equation for this type 
are discussed and implemented on a 'Forced' Lorenz system and a singularly perturbed system 
whose fast flow does not necessarily converge to an equilibrium. We explore two strategies. In 
one, we compute (locally) invariant manifolds of the fast dynamics, parameterized by the slow 
variables. In the other, the choice of our coarse variables automatically guarantees them to 
be 'slow' in a precise sense. This allows their evolution to be phrased in terms of averaging 
utilizing limit measures (probability distributions) of the fast flow. Coarse evolution equations 
are constructed based on these approaches and tested against coarse response of the 'micro- 
scopic' models. The second type of coarse variables are defined as (non-trivial) scalar state 
functions that are required by design to evolve autonomously, to the extent possible, with the 
goal of being candidate state functions for unambiguously initializable coarse dynamics. The 
question motivates a mathematical restatement in terms of a first-order PDE. A computational 
approximation is developed and tested on the Lorenz system and the Hald Hamiltonian system. 
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1. Introduction 

Obtaining coarse response of a system of ODEs containing rapidly oscillatory or decaying (stiff) 
response without detailed information on the evolution of the original (fine) variables is an 
interesting, but challenging task. The coarse response of the fine dynamics is sought with the 
following (ideal) objectives: (1) For oscillatory fine response, the fundamental frequency of 
the coarse dynamics should be orders of magnitude less than that of the fine dynamics. (2) 
The coarse dynamics should be autonomous in the sense that its evolution is independent of 
fine initial conditions. One possible application of our work would be in Molecular Dynamics 
(MD) simulations. MD is based on a particle description of atoms or molecules and treats the 
interaction between the particles using Newtonian dynamics. However, a key difficulty in MD is 
the extremely large separation between the time-scales of atomic bond vibrations (femtoseconds) 
and even the smallest time-scales of engineering interest (nano to micro seconds). This limits 
MD to unrealistically small time steps for understanding coarse evolution at engineering time- 
scales. In multi-scale analysis, it is beneficial if in the time integration of the system, the time 
step in the slow/coarse dynamics is much larger than in the atomistic dynamics. The objective 
of this work is to develop ideas towards the computational coarse-graining of such dynamical 
systems. In particular, we focus on methods for defining good 'slow' variables for probing 
macroscopic response and defining their closed evolution (i.e. in terms of themselves), to the 
extent possible. 

The first part of our work relates to defining the evolution of coarse variables. These objects 
are the 'slow' state variables in a system with an explicit split into slow and fast motions (e.g. 
constants of motion for the fast flow), and/or time averaged functions of fast variables. The 
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latter have the interesting property that their evolution is necessarily slow on the fast time scale, 
as shown in [3j. Thus, whether a-priori available, or induced by time-averaging, we deal with 
singularly perturbed systems of autonomous ODE. The approaches that we apply to define the 
coarse evolution equation are the following: 

• The first approach is to generate an augmented fine system with the introduction of time- 
delayed fine variables (where time delay is the averaging period) , and compute an appropriate 
collection of local/overflowing invariant manifolds of the discrete-delay augmented fine ODE 
system as graphs parameterized by the coarse variables a priori. With the computed collec- 
tion of invariant manifolds, a closed dynamics for the coarse variables can be defined. We 
refer to this procedure as Parameterized Locally Invariant Manifolds (PLIM) [U l22"j [3]. 
The invariant manifolds (computed locally) are considered as a representation of the fine 
variables on the coarse phase space. 

• The second question we explore originates from the reduced order approach (the Tikhonov 
approach [24J ) that utilizes the differential-algebraic equations (DAE) representing slow man- 
ifold^) to approximate the limit dynamics. The idea is to check whether, even when the 
fast flow does not converge to an equilibrium, the slow evolution of the first moment of the 
invariant measure generated by the fast flow for fixed slow variables is well approximated by 
the DAE of the Tikhonov method. 

• The third approach may be viewed as a natural extension of the method of averaging mo- 
tivated by practical considerations, naturally for systems whose limiting fast flows may not 
converge to an equilibrium. The theory of [8, 7j is applied here, where the limit behavior 
of the fast flow (parameterized by the slow variable) in a singulary perturbed system is 
characterized by an invariant measure rather than algebraic equations which are used in the 
Tikhonov approach. However, this theory requires that variables that can be averaged be 
'slow', i.e. orthogonal to the fast flow in a precise sense, but does not provide a prescription 
for generating such 'orthogonal observables'. Following |23j we develop and demonstrate a 
natural class of slow variables based on time averaged coarse observables, and generate their 
evolution, where this approach is referred as 'Practical Time Averaging'. 

The second part of this paper relates to the definition and computation of instantaneous 
(as opposed to time-averaged) coarse variables that aims to allow for an autonomous coarse 
response. Defining unambiguously initializable coarse dynamics with information only on the 
coarse state is the ideal goal (i.e. the time derivative of the coarse variable should be uniquely 
determined given a coarse state). However, this requirement for generating an autonomous 
coarse response can hardly be satisfied, as noted in the implementation of PLIM (for example, 
in the work of |22[ |4j) and the theory of [H [7J. In the case of PLIM, the invariant manifolds 
computed on the coarse phase space are not unique. A physical interpretation is that there could 
possibly exist multiple fine states that are associated with one single coarse state and these fine 
states affect the coarse evolution out of the attained coarse state. Thus, one requires knowledge 
of fine initial conditions to ensure a correct coarse response. In Practical Time Averaging, the 
'right-hand side' term of the coarse evolution equation is interpreted with an invariant measure 
which generally cannot be determined by the coarse state only. This second part of the paper 
develops further the preliminary ideas set forth in [2] in an attempt to define a class of coarse 
variables whose dynamical response can be uniquely defined by a given coarse state. 

This paper is organized as follows: In Section 2, we provide the theoretical framework of 
PLIM and present the numerical results from the implementation of PLIM on a 'Forced' Lorenz 
system. In Section 3, we show the procedures to write down the coarse evolution equation for 
singularly perturbed systems based on the Tikhonov approach and the Young measure theory 
(i.e. [HI [Tl [23]). Each of these approaches is followed by numerical studies of model problems. In 
Section 4, we formulate a first-order PDE to generate a class of coarse variables that formally 
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allow autonomous coarse response. We obtain such coarse variables from solving the PDE and 
discuss some desirable properties of these coarse variables. We end in Section 5 with some 
discussions of the strengths and drawbacks of these multi-time scale modeling techniques. 



2. Parameterized Locally Invariant Manifolds 



The main ingredients of PLIM are discussed in [TJ HJ l22"j |3] . PLIM is based on more-or-less stan- 
dard ideas for computing invariant manifolds using the invariance equation of an autonomous 
ODE system |21[ \§\ \TU[ \TE[ fTB"l [H] with one essential, and important, difference. The coarse 
'retained' master variables whose dynamics is sought are not some of the fine variables, or a 
function of those. Instead, they are running time averages of functions defined on fine trajec- 
tories. With a simple augmentation of the fine dynamics as proposed in [3], a coarse dynamics 
can be obtained for the time averages. The advantage of this strategy lies in the fact that 
time-averaging provides a natural way to separate time-scales in the system consisting of the 
augmented fine ODEs and the coarse variables, even though the original (non-augmented) fine 
ODE system may not come equipped with an in-built and identifiable separation of scales in its 
variables [3]. 



2.1 THEORETICAL FRAMEWORK 



We are interested in nonlinear, autonomous fine dynamics defined as 

ft (<> = *(/(<).'(»> ; |w = f £ ««» (i) 

/(0) = /» ; /(0) = J„ 

where / is an A-dimensional vector of fine degrees of freedom and H is a generally nonlinear 
function of fine states, referred to as the vector field of the fine dynamical system. N can be 
large in principle, and the function H rapidly oscillating. / is an n dimensional vector of slow 
loading variables that do not vary appreciably compared with the oscillation of fine variables. 
Here, T is a non-dimensional parameter that represents the ratio of natural frequency of the fine 
dynamics to the loading frequency, i.e. T := w//cj, where cjf denotes a characteristic frequency 
of the fast dynamics and uj is the time scale of the loading (e.g. constant rate of monotonic 
loading, frequency of cyclic loading). L takes the form L(I) := L(I)u)f, where L is a bounded 
smooth function with dimensions of the load vector /. For Uf fixed (which is what we have in 
mind), by definition, T — > oo as u — > 0. 

Let c denote the coarse variables of dimension m and A be a user-specified function of the 
fine states producing vectors with m components whose time averages over intervals of period 
r can be measured in principle and are of physical interest. Given the fixed time interval r 
characterizing the resolution of coarse measurements in time, a coarse trajectory corresponding 
to each fine trajectory /(•) is defined as the following running time average: 

c(t):=- f +T A(f(p))dp, (2) 

T Jt 

where r < 1/uj. Thus, r is the period of time averaging, referred to the fast time scale. Clearly, 



d ^it)= l -[K{f{t + T))-K{f{t))\- (3) 
From the definition of time-averaging, forward trajectories //(•) and !/(•) are introduced as 
ff(t):=f(t + r) 

I f (t):=I(t + r), {) 
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then 

d ll(t) = H(f f (t),I f (t)) 

§<*) = **(/,(*)). 
With these definitions in hand, an augmented system is introduced: 

d ll(t)=H(f f (t),I f (t)) ; f f (0) = U 

f t (t)=H(f(t),I(t)) ; /(0) = /, 

f(t) = ±L(I f (t)) ; 7,(0) = /« (6) 

^(t) = ^(/(t)) ; 1(0) = /, 
oV 1 

-(t) = -[A(/ / (t))-A(/(t))] ; c(0) = c*, 

where we assume that /* and /** are arbitrarily specifiable. Thus, the set of evolutions described 
by ^ is strictly larger than that defined by ([!]), the latter being a subset of the former. 

Now we seek functions G^(c, I, If) and G l (c, I, If), i = 1, . . . , N that satisfy the first order, 
quasilinear partial differential equations 

™dG J f ( l , , \ "lfdGi . dGj . \ 

Efl/ - A*(G)] j + £- ( L*(I) + ^(It)) = H (Gf, If) 



fc=i j=i \ j 




(7) 



at least locally in c-space. These are generated by assuming Gf = ff and G = f, and substi- 
tuting in ([6]). Suppose we have such a pair of functions over the coarse domain containing the 
point 

c* := c(0) (8) 
defined from ([!]) and Q which, moreover, satisfies the conditions 

G f(c*) = /** ^ 

Also, given an initial state (/*, /*) we denote (/**, I**) as the solution of ([l]) evaluated at time 
r. It is easy to see that a local- in-time fine trajectory defined by 

Tf(t) := G f (c(t)) 

T(t) := G(c(t)) 1 ; 

through the coarse local trajectory satisfying 
dc 1 

- = -[A (G f (c, I, If))- A (G(c, J, J,))] (n) 
+ the evolutionary equations for 7, 2j 



is the solution of ([6]) (locally). A solution pair (Gf,G) of ([7j) represents a parametrization of a 
locally invariant manifold of the dynamics ([6]). Thus, ( |11[ ) is the closed coarse evolution corre- 
sponding to the augmented fine system ^ and ([9]), £mi on/?/ locally. 
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2.2 A MODEL PROBLEM: 'MONOTONICALLY FORCED' LORENZ SYSTEM 

We test the idea of PLIM on the Lorenz system |17j , and take the time average of state variables 
in the fine phase space as coarse variables. As is known, the Lorenz system is a three-dimensional 
chaotic system evolving over time in a complex, non-repeating pattern within a bounded region 
of phase space (for appropriate parameter values). The ODE is given as: 

x = a(y — x) 
y = 7X — y — xz 

a (12) 

z = xy — pz 

a = 10, p = 8/3, 7 = 25. 

We are interested in the evolution of running time averages of x over different periods r, while 
the Lorenz system contains high frequency non-periodic bounded oscillations. These oscillations 
more or less average out to zero. Thus, for the purpose of testing our ideas, we construct a 
'Forced' Lorenz system whose fine variables show high frequency oscillations superposed on a 
monotonically increasing/decreasing profile. 

Let x(t) denote the trajectory whose time average has a monotonically increasing feature 
while keeping the high frequency oscillations at the fine scale. Let L(t) be a trajectory that is 
monotonically increasing. One possible construction of x would be: 

x(t) = x{t) + L(t), where L = ^——r. (13) 

1 + L 

The trajectory with a general increasing pattern can be obtained by simply considering a linear 
combination of the original trajectory with a monotonically increasing curve. However, to derive 



the corresponding ODE requires some more work. From (13), we have 

x(t) = x{t)-L{t). (14) 
We also have 

x = x + L = a(y — x) + -. (15) 

1 + L 



Substituting (14) into (15) and (12)2,3, we obtain the ODE of the 'Forced' Lorenz system: 

1 



x = a(y — x) + aL + 

j. -r u 

y = 7X — y — xz — 7L + Lz 

- R T (16) 

z = xy — pz — Ly 

1 

L 



l + L 

For simplicity, we denote x as x and rewrite the ODE as the 'Forced' Lorenz system: 

1 



x = a{y — x) + aL + 

j. -r u 

y = 7 x -y - xz - 7L + Lz 

z = xy — j3z — Ly 
1 

L 



l + L 

Figure [T] shows the trajectories of the original Lorenz system and the 'Forced' Lorenz system in 
the (x, y, z) phase space. The initial conditions are set to be the same. For the Lorenz system, 
the trajectories evolve around the steady state points with no preference of moving towards any 
direction, while for the 'Forced' Lorenz system, the trajectories move towards the positive x 
direction over time, keeping the same amplitude in the y and z directions. 
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We take the running time average of x as a coarse variable: 



x(t) 



t+T 



x(p)dp =>• x{t) = - [x(t + r) — x(t)] 



i x f 



where Xf:=x(t + r). (18) 



We want r to be large enough. The highly oscillatory curves can be reduced to a constant 
by taking an infinitely long time average. Practically, the time averaged variable will be grad- 
ually smoothed by increasing the time average, as shown in Figure [2] The estimated lowest 
fundamental period of the Lorenz system is around 0.0455 (This approximation is based on 



calculating the eigenvalues of the linearized system of (12). Note that time is dimensionless in 



the original definition of Lorenz system |17j . so this estimated number is also dimensionless). 
The time step for evolving the fine dynamics is of the order of 10~ 3 . We set the averaging 
period to be 50, which is 5 * 10 4 times larger than the fine time step. 




x 

■x(f = 10) 
-x( T = 30) 
■x(f = 50) 



Figure 2: Comparison of x(t) and x(t) with different time averages 



With the introduction of the time shifted variables Xf, yj, zj and Lf, the corresponding aug- 
mented system takes the form 

x f = a(y f - x f ) + aL f 



1 + Lf 

y'f = i x f -vf- x / z f - i L f + L f z f 

z 'f = x fVf - /% - L fVf 



(19) 



L 



f 



1 

l + L 



f 
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Taking x and L as the coarse variables, we first construct a mapping on the coarse domain 
that takes values in the fine phase space, i.e. x = Gi(x,L), y = G2(x,L), z = Gs(x,L), 
Xf = G±(x, L), yf = G§(x, L), zj = Gq(x, L) and Lf = Gt(x, L). Thus, the governing equations 
for solving the invariant manifolds have the form 

<>C ' " (T ' J N + ^ (t^Tt) = Hi(G) V i from 1 to 7, (20) 



dx V t J dL \1 + L, 
where the right-hand side terms are 

Hi(G) = a(G i+l - d) + aG 7+{i _ 1)/3 + 1 ; i = 1 or 4 

~ 1 + (j 7+(i-l)/3 

Hi{G) = 7Gj_i - Gi — Gi-iGi + i — r fG 7+ / i _2)/3 + G f i+iG 7+ (j_ 2 )/3; i = 2or 5 

(21) 

iJj(G) = Gi-. 2 Gi-i - (3d - Gi_iG 7+ (j_3)/3; i = 3 or 6 

The coarse theory is defined as 

j_ G 4 (x,L)-G 1 (x,L) 
x = 

1 T ( 22 ) 

L = TTI- 

To obtain the time average of x in the original Lorenz system, we simply subtract the function 



L from x computed from ( 22 ) 



Figure [3^a) shows the comparison between x(t) from evolving the fine theory and the coarse 
theory with different coarse-to-fine time-step ratios (denoted by c/f). The coarse trajectories 
exhibit consistency of a general growing trend, although the two curves do not perfectly agree. 
When the coarse responses with different time-step ratios are compared, the results agree well 
up to c/f = 10 4 . This observation is consistent with the previous choice of time average, which 
is of the order of 10 4 of the fine scale. Figure [3^b) is a magnified version of the results from 
t=0 to 10, which demonstrates that the responses are much more consistent at the early steps 
of evolution. Figure [3](c) shows a collection of coarse trajectories computed by evolving the 
fine theory, where the difference of any two coarse initial conditions defined in the Z 2 -norm is 
smaller than 10~ 3 , while the difference of the associated fine initial conditions could be quite 
large. At the coarse scale, the differences among the coarse initial values are considered trivial, 
nevertheless the trajectories do not totally overlap due to the chaotic nature of the Lorenz 
system. On the other hand, all of the trajectories follow a gradual increasing pattern as the 
way it is built up for the 'Forced' Lorenz system. The coarse responses from evolving the coarse 
theory show a better consistency, as seen in Figure [3^d), where the trajectories are almost 
overlapping. Here the initial conditions are the same as before. This phenomenon implies 
coarse stability of the constructed coarse model via time-averaging in that small differences of 
coarse initial values lead to small differences in coarse evolution. The chaotic feature observed 
at the fine scale is no longer maintained at the coarse scale. 

Figure [4] shows a group of coarse trajectories from evolving the fine theory, where the initial 
data are arbitrarily chosen, that is to say, the fine initial conditions corresponding to each coarse 
trajectory can be significantly different. However, due to the long time averaging, the coarse 
responses exhibit a similar increasing pattern as noticed in previous results. 
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Figure 3: (a) Comparison of coarse and exact response for different coarse-to-fme time-step 
ratios; (b) Magnified version of (a) from t=0 to 10; (c) x(t) from fine theory with different fine 
ICs, where the Z 2 -norm of any two corresponding coarse ICs is within 10 -3 ; (d) x(t) from coarse 
theory, where the coarse ICs are the same as (c) 




Figure 4: x{t) from fine theory with arbitrary fine ICs 



In conclusion, we succeed in obtaining consistent coarse response from solving the coarse theory 
with much larger time-step. We learn that by taking a large enough averaging period, we can 
obtain a coarse dynamics that is insensitive to fine ICs, while retrieving some key features at 
the coarse level; on the other hand, the coarse dynamics is not aimed to capture the chaotic 
feature and high frequency dynamics at the fine scale. With small differences of coarse ICs, the 
coarse responses show no difference. Both of these properties are desirable since the dynamics 
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at much larger time scale than that of fine dynamics is of interest to us. 
3. Approximation of High Frequency Oscillatory Response 

The problems of interest here are singularly perturbed systems that display highly oscillatory 
response and whose limiting behavior often retains fast oscillatory response. Many physical 
problems can be expected to fall in this category. For instance, in MD simulations where exter- 
nal loading is applied to the system, the loading rate is usually at least 10 9 times slower than 
the fundamental frequency of atomic vibrations. This large separation of time-scales leads to 
singular perturbation in the fine system ([!]) since 1/T becomes a small number. In this study, 
we are interested in understanding the evolution of time averages of phase functions of the 
fast flow on the slow time scale of loading. Two approaches are developed to approximate the 
time averaged dynamics. The first approach utilizes the Tikhonov method |24| of differential- 
algebraic equations (DAE). The second approach is a natural, practical extension of the method 
of averaging proposed in [8j [7| . We construct the evolutionary equations of time averaged quan- 
tities following the work of [23] and refer to this collection of theoretical and numerical ideas as 
'Practical Time Averaging'. 

3.1 THE TIKHONOV APPROACH 

The system we are concerned with has the same form 

df 
dt 

dl 1 



dt T i(/ » < 23 > 



1 



t+T 



c(t) = - / A(/(p))rfp, 

T Jt 

where / G R N , I € M n and c E R m . Recall that if ijj is the loading frequency, then T — > oo 
as oj — > 0. Define a new time scale s := t/T and denote e := 1/T. Note that since T is 
dimensionless, s has the unit of time and e is a dimensionless number. For the functions / and 
I of the fast time variable t, we define the functions / and I of the slow time variable s as 
follows: 

f(s) = f(sT), I(s) = I(sT). (24) 
Then, on the slow time scale the system reads 

ef = H(f,L) 

ds (25) 

The coarse variable is defined as 

c(s) = c(sT) = - f T A(f(p))dp, (26) 

T JsT 

and we assume that r = n/u with n a non-dimensional constant and < k < 1. Introducing 
the change of variables r = p/T, we have 



1 



s+X 



c(s) = - I A(/(r))dr, (27) 
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where A = t/T. Therefore, A = k/(ujT) = n/cjf is a constant as e — > 0, with physical dimensions 
of time. 

Removing all overhead tildes, we have the slow-fast system, posed on the slow time scale: 



4, 

dl 



ds 



L{I) 



(28) 



l r s+x 

c(s) = -J A(f(r))dr. 



Let / be the solution of the following differential-algebraic system obtained from (28) when e 
is set to be equal to 0, namely 



= H(f,I) 
dl 



ds 



L(I). 



(29) 



Then by using the chain rule, we obtain 



dH df 
df ds 



dH dl 
~dTds~' 



(30) 



So the new evolution equation (rate-independent) for the fine variables is 



df 
ds 



(dH 



dH dl 



dl ds 

The new system becomes 

df _ (dH\~ l dH 
~d~s~ ~\dj) ~dT 



(31) 



ds 
c(s) 



(32) 



s+X 



A(/(r))dr 



with initial conditions on / satisfying H(f, I) = 0. 



A classical result of Tikhonov [24J states the solutions of (28) converge, as e tends to zero, to 



the solution of the differential- algebraic system (29) under appropriate assumptions. A crucial 



condition for this result to hold is that for each fixed /, solutions of the differential equation 
df jdt = H(f, I) should converge, as t — >■ oo, to a stationary point /(I) such that H(f(I), I) = 0. 
However, a more general theory given by [H [5] shows that this condition may not be always 
satisfied. Instead, the convergence to the limit is in a weaker sense (in the sense of Young 
measures), namely, one gets information only about the probability measure of the fast flow. 
Our work is to test whether the slow time evolution of the first moment of the probability mea- 
sure of the fast flow parametrized by I(s) corresponds in any sense to the solution of the DAE 



of the Tikhonov approach. Realizing that (29) and (32) are equivalent given the ODE starts 



from an equilibrium point of the fast dynamics, we refer to the ODE system (32) as the DAE 



dynamics of the slow- fast system (28). The system, in general, is not well-defined when dH /df 
is singular. However, well-posed evolution can exist in many situations. At the practical level, 
this may be solved by trying different numerical algorithms (for example, using implicit Euler 
backward method on the points near the singularity), but we avoid an extensive consideration 
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of this fact in the following model problems. 



3.1.1 Singularly Perturbed System with Measure- valued Limit 

We discuss a specific realization of Example 5.1 in [5j. The slow-fast system is given as 

yi 

£t = V2 
as 

2/2 , \ 

£t = 9{x,yi,y 2 ) 
as 



(33) 



dx 
ds 



2/1, 



with state variables y%, y 2 and x. The function g (a;, 2/1,2/2) is defined as 

g(x,yi,y 2 ) = a(x){l - A(x)(yi - z up (x)) 2 )y 2 - {yi - z up (x)), 
where we set a(x) = (x — 3~22) 2 , A(x) = l/(x — 3~22) and 



x) = : -x + yi - yf = V x G K and 2/1 > 3 a }. 



(34) 



(35) 



for the dynamics at the upper branch (see Figure [5]). For the lower branch, we substitute 
a(x) = (x — (— 3~ 2 2)) 2 , A(x) = l/(x— (— 3~5 2)) and replace z up (x) with zi a (x) which is defined 
as 



z lo(x) = {yi : -x + yi - yf = V x £ E and y x < -3 2 }. 



(36) 



Thus the corresponding differential-algebraic system for (33) has the form 



= 2/2 

= -x + yi - yf 
dx 
ds 



(37) 



2/1- 



However, for e — >■ the limit behavior of (33) does not follow the DAE (37). Instead, for each 



fixed x, the fast motion will converge to periodic orbits (also referred to as limit cycles) around 
either the upper or lower branch of the equilibria profile given by ( 37 ) , as shown in Figure [5j 
The blue curve exhibits the trajectory from running the slow- fast system (33), plotted on x-yi 
phase space. Meanwhile, there are fast motions around y 2 = (not shown) during evolution. 
The equilibria profile is shown in red with y 2 equal to zero constantly. We start the dynamics 
at x = —1.875. The initial values of (2/1,2/2) are set to be (1.5,0). However, in fact, starting 
from any values of (2/1,2/2) with 2/1 > z up {x) will lead the fast motions to converge to the limit 
cycle around the equilibrium point within a short time interval with x hardly changing. The 
trajectory continues its fast motion in 2/1-2/2 plane close to the limit cycles (even in the limit 
e — > 0), never reaching the slow manifold defined by g(x,y\,y 2 ) = with a slow movement 

3 

along the positive x direction till x = 32 2. Similarly, when starting from a point on the lower 
branch, say x = 1.875, the fast solution will converge to the limit cycle shortly following a slow 
drift along the negative x direction till the point x = —32 2. 



The DAE dynamics for (33) according to (32) is 



dyi 
ds 
dy 2 
ds 
dx 

ds 



dg_ 
dx 



dg 
dyi 



2/1 



(38) 



2/1- 



12 



The numerical result is shown in Figure |6j Initial conditions are chosen to satisfy H = 0. There 
are two singularity points: (x, 2/1,2/2) = (32 2,3~2 ; 0) and (21,2/1,2/2) = (— 3s 2, — 3~ 2 , 0). We 
start the simulation from the point ( — 1.875,1.5,0) on the upper branch until the trajectory 

3 1 

reaches the singularity point (32"2,3 _ 2 ; 0), and then repeat the same procedure on the lower 

3 1 

branch till the point (—32 2, — 3~2 ; 0). 2/2 is constant zero in both cases. It can be seen that the 
time averaged behavior matches well with the DAE dynamics. The two trajectories are almost 
overlapping. Although the DAE is not the appropriate limit solutions of the system, it succeeds 
in capturing the time averaged behavior of the limit dynamics. 




-2 -1J -1 -0.5 0.5 1 1.5 2 



Figure 5: Diagram of equilibria and limit cycles. Red curve indicates the equilibria of the fast system 
and blue curve shows the limit cycles generated from the slow-fast system in the limit e — > 0. 
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X 

Figure 6: Comparison of time averaged dynamics and DAE dynamics. Blue curve indicates y 1 vs. x; 
red curve indicates y% vs. x. Solid line indicates stable dynamics; dash line indicates fast motion of limit 
dynamics. 



3.1.2 'Oscillatory Forced' Lorenz System 



The original Lorenz system is a three-dimensional chaotic system. The trajectories tend to 
converge on a topologically complicated, but bounded set that can be enveloped within another 
bounded and contiguous region of 3D phase space. We apply the same strategy as in Section 
2.2 to make the averaged x variable evolve sinusoidally with loading period while keeping the 
high frequency oscillation. By doing so, we assure that the trajectories of the 'Forced' Lorenz 
system will remain bounded (under the condition that the loading is bounded) and exhibit 
chaotic patterns along a slow sinusoid. Meanwhile, the 'Forced' Lorenz system will show an 
explicit separation of slow and fast motions. Towards this goal, we define 

x(t) = x(t) + L(t) 

w : ' w (39) 

where L = ojg and g = —uiL. 
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Obviously uj indicates the loading frequency. Prom (39), we have 

x(t) = x(t) - L(t) 
and 

x = x + L. 



(40) 



(41) 



(42) 



Substitute (39) and (41) into the original Lorenz system (12), and remove the overhead tilde 
notation, 

x = a(y — x) + oL + ujg 

y = 7X — y — xz — 7L + Lz 

z = xy — f3z — Ly 

L = ug 

g = -ujL. 

Assume s = ujt, the system on the slow time scale takes the form 
dx 

lo— = a(y — x) + aL + ujg 
as 

dy 

uj— = jx — y — xz — jL + Lz 
ds 

uj— = xy - f3z - Ly 
ds 

dL 



(43) 



ds 
dg 

ds 



-L. 



For uj — > 0, the equilibrium points of the fast system at s = are denoted as EP1 
EP2 = (-8, -8, 24) and EP3 = (0, 0, 0). 
The DAE dynamics following d29|)-(|32|) is 



5,24), 



dx 

ds 

dy 

ds 

dz 

ds 

dL 

ds 
dg 

ds 



ujL[(L — x)(x — L) — P] 



a(-2Lx + L 2 + x 1 - /?7 + f3z + f3 - Ly + xy) 
ujL(Ly - (3z + ^7 - xy) 
a(-2Lx + L 2 + x 2 - /3j + fiz + /3 - Ly + xy) 

ujL(y + 7X — 7L — zx + zL) 
a(-2Lx + L 2 + x 2 - £7 + fiz + (3 - Ly + xy) 



-L. 



+ 9 



(44) 



Since uj exists in the right hand side terms in (43), it remains in the DAE dynamics. With 
further assumption that uj — > in (44), we obtain a simple form of the DAE dynamics: 

dx 
ds 
dy 

ds 

dz 

ds 

dL _ 
ds ^ 



9 





(45) 
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Figure [7] shows the fine variables and their time averages from evolving the system ( 43 ) . The 
evolution of x is a superposition of the high frequency oscillation with the periodic loading 
whose time average coincides with the sinusoidal loading. The time averages of y and z are 
nearly constants, despite the oscillatory motion on the fast time scale. Figure |8^a),[8Fb) and[8^c) 
display the time averages of x, y and z, respectively, starting from different initial conditions. 
It is noted that the averaged behavior, or the dynamics on the slow time scale, is insensitive to 
the fine initial conditions (for the cases tested). The trajectories that start from different points 
converge to a unique time averaged dynamics (this statement is limited to the numerical tests 
we did) and follow the time averaged dynamics consistently. 

However, the DAE dynamics (45) does not give consistent results. The trajectories starting 
from equilibria do not coincide with the time averaged dynamics, as seen in Figure |9j Fig- 
ure[9^a), (b) and (c) present the results from fine dynamics and DAE dynamics of the variables 
x, y and z, respectively. The black curve denotes the time average from the fine theory. Al- 
though the trajectories evolve with the same pattern, they do not agree with the time averaged 
behavior since the initial conditions corresponding to equilibria of fast dynamics do not coincide 
with the averaged dynamics. The last test is to set initial conditions for the DAE dynamics as 
the average of (x,y,z) from previously assigned fine initial conditions. In this case, (x, y, z) are 
consistent with (x, y , z) (see the magenta curve in Figure [9]) . 

There is an analogue between the 'Forced' Lorenz dynamics and Example 5.1 discussed in 
Section 3.1.1 as can be seen from a plot like Figure [7]^a). The fast motion occurs in (x,y,z) 
space (the dynamics in the y and z directions is not shown) with x oscillating around the slow 
sinusoid and (y, z) moving around a fixed point close to (0, 20.6). 
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Figure 7: Comparison of fine evolution and corresponding time averaged dynamics 
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Figure 8: Time averaged dynamics starting from various initial conditions 
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Figure 9: Comparison of time averaged dynamics from fine and DAE dynamics starting from 
initial conditions 
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3.2 PRACTICAL TIME AVERAGING 



Different from the Tikhonov approach, but including it when appropriate, the work of [7] 
building on earlier work of [8] allows dealing with limits of fast flows that are not necessarily 
equilibria. The theory accounts for the limit behavior of fast flows that may be supported on 
w-limit sets of complicated topology. We state the relevant part of the theorem below, following 
the terminology presented in [23]: 



Suppose the dynamical system of interest, written on the slow time-scale, is (28). Let f e ^'f° 
represent solutions of the fast system 



df 



E,I,fQ 



ds 



e,I,fo 



I) ; f e ' IJo (o) = f 



(46) 



with / fixed. Rewrite the fast system with the introduction of the fast time scale t = s/e, 

d f(0),I,f 



dt 



-(t)=H(fW> I 'f°(t),I); /(°)^/o (0) = /o 



(47) 



where 



/ 



(0),/,/o 



(*) = / 



(0),/,/o 



/ 



s . 



(48) 



Assume that for < t < oo, f 6 ' 1 ^ lie in a bounded set for each /. Hence, on [0, T] we will have 
for any continuous function F 



F(f' IJ »(s)) -> / F(7K J)/o ( 7 )d7 weak * L°°[0,T] 



(49) 



where (J> s I f is the Young measure generated by the sequence f 6 ' 1 ^ £ L°°(0, T; M. N ). Moreover, 
the theory proves that /U s i/ is an invariant measure supported on the w-limit set of /(°)> J >/°. 
Note that f 6 ' 1 ^ is obtained from running the fast system (47) for a sufficiently long period of 



time, starting from /q while keeping / fixed; The function F in (49) can be linear or nonlinear 



In situation when the w-limit set of the fast dynamics is not an equilibrium point, the Tikhonov 
approach may not in general be expected to approximate the time evolution of the first-moment 
of nonlinear functions of state. Furthermore, the theory proves that coarse evolution equations 
can be written down for state variables that are 'slow', i.e. orthogonal to the fast flow in a 
precise sense, but does not provide a prescription for generating such 'orthogonal observables'. 

We follow the work of [23], which provides a straightforward way of generating a large 
class of slow variables via time averaging that are shown to be 'orthogonal observables' (for a 
corresponding infinite dimensional dynamics) for which a slow dynamics can be written down. 
We refer to this approach and the numerical approximation ideas we develop in this paper as 
'Practical Time Averaging'. 



Recall the coarse variable defined on the slow time scale (28) and write down its evolution: 



c{s) 



1 



s+X 



A(/(r))dr 



dr 1 

-(s) = -[A(f( S + X))-A(f(s))}. 



Let us define the sequence (on e) of smooth functions c £, ^° such that 



1 



s+A 



c £ >f»(s) = ± / A(f £ ' I( - r ^f°(r))dr, 



(50) 



(51) 



where we follow the notation in ( 46 ) but now / is a function of the slow time. Then the evolution 
equation is 

dc £ ' fo 



<h {s) = \ [A(f £ ' Hs+x)Jo (s + A)) - Aif^^is)) 



(52) 
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Let c*° denote the weak limit of the sequence c e, f°, and using (49), we have 

dc? _ 1 
ds A 



A(7)M; 



+X,l(s+X),f 



A(7)^,/( s ),/o(7)^7 



(53) 



in the weak (variational) sense. (53) can be considered as the instantaneous evolution equation 
of the limit slow variable of the sequence c e, ^°. It should be noted that if any arbitrary 
instantaneous function, say A(f), is taken as the coarse variable, we have 



-VfA-H. 



(54) 



ds £ 

Since 1/e is built in the evolution equation, the Young measure limit is usually not applicable 
unless the function A(f) satisfies VfA-H = 0. However, our definition of coarse variables allows 
constructing coarse evolution equation for more general kinds of averaged phase functions. For 
example, to approximate the change of temperature of a MD system, where temperature can 
be interpreted as a function of particle velocities but normally does not satisfy V/A • H = 0, we 
can define the time average of temperature and follow the procedure from ( 50 ) to ( 53 ) to write 
down the coarse evolution equation. 



To numerically approximate the coarse dynamics ( 53 ) , we can apply a three-stage procedure 
suggested in [BJ: 



1. Evolve the fast system (47) with slow variables fixed at time step s and s + A, respectively 



to generate the Young measures as averages of M Dirac masses at M values of /, i.e. 

1 M 

A(7K/ {s) ,/o(7)d7 « m E A(/ (0) ' /(s) ' /0 ) 



1 M 

/ A( 7 K w(s+A)i/0 (7)d7 « T. E A (/ (0) * 7 



(55) 



(s+A),/ 



where M is chosen to be large enough for the averages to converge. 

2. Extrapolate to obtain the coarse variable at the next coarse step s + T using the Forward 
Euler method, 



c fo (s + T)^c }0 {s) + 



fo l 



A 



A(7K+A,/(s+A),/ (7)d7 



A(7)/V(a),/b(7)<*7 



(56) 



where T is the coarse time step. Usually T is set to be much larger than A. I(s + T) can be 



obtained from evolving (28)2 solely. 



3. Reconstruct to find fine data that matches c^°(s + T) and I(s + T). 

4. Repeat. 

The primary shortcoming of all of this theory as acknowledged in [HI [7] and carries over for the 
model with time averages in |23| is that the theory is silent about fine initial conditions required 
to generate the Young measures at the discrete slow time- levels. This is a serious practical issue 
and renders fo as an essentially free parameter whose choice, nevertheless, affects results and 
has to be made on an ad-hoc basis based on numerical experimentation. We mention here that 
the unsubstantiated assumption of ergodicity is often made to sweep this important issue under 
the rug. 

Given the above uncertainty and the need to use it twice to follow the evolution prescribed 



by (56), an alternative procedure is the following. Based on (48), the sequence of functions 



c £, ^°(s) can be written as 



00 



,(0),/o 



.(0),/o 



(*)• 



(57) 
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Using the substitutions p = r/e, r 
interpreted on the fast time scale as 



A/e and (48), the right-hand-side term in (51) can be 



s+\ 



r))dr 



t + T 



A(/(°>'W°(p))dp, 



(58) 



where we make the realistic assumption that the loads I do not vary appreciably in the (fast- 



time) interval [t,t + r]. From (57) and (58), we have 



„(0),/o 



t+T 



A(/(o)^(t),/o (p))dp) 



(59) 



which is essentially of course. We are interested in c £ 'f° as e — > 0, which by definition, is the 
time average of functions of the fast variables as r — > oo when the slow variable I is fixed. We 



simply follow (59) to calculate this average at each coarse time step T. Theoretically, we need 



to run the fast system for an infinite period of time. However, for the numerical tests conducted 
here we have seen that the average converges quickly, thus taking a reasonable period of time 
is sufficient. 



Figures 10 and 11 present the time averaged variables from running the slow-fast system and 
the limit slow variables, of the model problems in Section 3.1.1 and 3.1.2, respectively. For both 
cases, we observe good agreement between the time averaged dynamics and the averaged limit 
dynamics. It is worth mentioning that: 

1. For the 'Forced' Lorenz system, we successfully approximate the coarse dynamics starting 
from 'arbitrary' fine initial conditions. The coarse dynamics is insensitive to the fine initial 
conditions, namely, starting from a point that is far away from the w-limit set, the fast 
dynamics will jump to the set quickly and evolve around it. 

2. We also test x 2 as a nonlinear function of the fast variables, where the result is consistent 
with x 2 , as shown in Figure 11 d); In contrast, using the Tikhonov approach and defining 

x 1 



? 2 as a coarse variable of interest, we will finally obtain the limit dynamics of x 2 



rather than 



x- 



3. As e — > 0, the fast dynamics of the problem (33) converges to periodic limit cycles, and 



the fast dynamics of the 'Forced' Lorenz system exhibits chaotic patterns; However, in both 
cases, the time averaged coarse variables approximate the time average of the limit solutions 
(which retain fast motions). 



4. For the problem in Section 3.1.1, the limit cycles in the upper and lower branch center around 
z up (x) and zi (x), respectively, where the centers satisfy H = and the time averaged 
trajectory and the DAE dynamics follows this slow manifold of the line of centers; the 
averaged trajectory of the 'Forced' Lorenz system is distinct for the three initial conditions 
that satisfy fast equilibrium, which causes the mismatch between the time averaged dynamics 
and the DAE dynamics starting from H = 0. 
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Figure 10: Comparison of time averaged dynamics with PTA dynamics. Blue curve indicates y 1 vs. x; 
red curve indicates yi vs. x. Solid line indicates slow dynamics; dash line indicates fast dynamics. 
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Figure 11: Comparison of time averaged (coarse) variables from fine and from Practical Time Averaging 
(PTA) 



To show the necessity of defining 'orthogonal observables' as coarse variables for the Young 
measure theory to be applicable, we follow (54) to write down the coarse evolution equations 
for instantaneous functions that are generally not orthogonal. We test this idea on 'Forced' 
Lorenz system with the chosen coarse variables being x, y, z and x 1 . The corresponding coarse 
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evolution equations are 
dx 1 



ds 
dy 



dz 
ds uj 
dx 2 



(a(y - x) + uL) + g 
(jx — y — xz — 7L + Lz) 
(xy - /3z - Ly) 



(60) 



1 



2x[-{a(y-x) + aL)+g 
ds \u 



dL 

ds 



9 



dg 
ds 



-L. 



The results from running the system of ODEs (60) are shown in Figure 12 The instantaneous 
functions x, y, z and x 2 are quite oscillatory due to 1/uj in the evolution equation. This high 
oscillation is inappropriate for a physical coarse variable whose evolution should be amenable to 
easy approximation on the slow time scale. On the other hand, the time averaged variables x, 
y, z and x 2 are much smoother and seem to be reasonable slow variables. This is one example 
of how PTA has the capability of producing robust slow variables for coarse dynamics. 
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Figure 12: Comparison of non-orthogonal instantaneous variables and time averaged variables on the 
slow time scale 



4. Variables for autonomous coarse response 

In the construction of coarse theory for fine systems, autonomous coarse response is usually 
desirable. However, the coarse theory based on the idea of PLIM or Practical Time Averaging 
cannot be guaranteed to be autonomous, i.e., the right-hand side term in the coarse theory 
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is not uniquely denned from one single coarse state. The objective of this part is to explore 
this question and mathematically formulate a PDE for a class of coarse variables to formally 
allow an autonomous coarse response, where we further explore some appealing features of the 
selected coarse variables. 

4.1 MATHEMATICAL SCHEME 

Suppose the problem of interest is 

f = (61) 
/(0) = /„ 

where / £ M. N . Based on the work in [2], we define the coarse variable c(t) as an instantaneous 
scalar function on the fine phase space and denote it as n. Thus we have 

dn 



it) := n(/(t)) => c(t) 



Of 



(f(t))H(f(t)). 



(62) 



From previous work [22] on constructing the coarse dynamics of the Lorenz system using PLIM, 
we realize that for an arbitrary coarse function n, the coarse dynamics cannot be ensured to be 



autonomous. Figure 13 shows the fine-to-coarse conversion of the Lorenz system by arbitrarily 
choosing x and z as coarse variables. Figure [l3^a) displays the fine trajectory on (x, y, z) space. 
The coarse trajectory is simply the projection of the fine trajectory on x-z plane, where we 
observe self-intersections of the trajectory, which simply imply a history-dependent behavior of 
the coarse evolution. 
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Figure 13: Fine to coarse conversion (a) Fine response; (b) Coarse response 



Therefore, the coarse variable needs to be specifically designed to formally ensure an autonomous 
coarse response. As outlined in [2j, a coarse variable is sought such that, if two fine states 
correspond to an identical coarse value, so does the coarse evolution from these fine states. 
Mathematically, this requires the existence of a function n for which there is a non-empty set 
of fine states defined by ( p3| ), where c is a constant and A c is another constant dependent on c 
(and independent of fine states), i.e. 

{/:H(/)=c and ^(f)H(f) = A c }. (63) 

The existence of such a function (with enough smoothness) would imply autonomous response 



by (62). Geometrically (63) implies that a level set of n(-) is also a level set of dH/df (■)!!(■). 
Therefore 



(64) 
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where A is an arbitrary scalar. By taking A as a derivative of an arbitrary function <3?(n) of the 
coarse variable, we have 



A := 



<9$(n) 
an 



d_ 
df 



0. 



One sufficient condition of (65) is to make II satisfy the first order PDE 

•TO -fir. 



(65) 



(66) 



Thus, a function II satisfying (66) satisfies the necessary condition for autonomous response of 
II, i.e. (64). In order to test this idea, we choose <E> to be an affine function of II, namely, 

dU 



df 



H = a + bU, 



(67) 



where a and b are constants. (67) is solved by implementing the finite element scheme. First, 
the coarse variable is discretized on the fine phase space, which takes the form 

N 



n(/) = DnV(/). 



(68) 



A=l 



The residual of (67) is minimized, which leads to the following linear equations 

N 



k AB * U B = b A , 



(69) 



B=l 
where 

kAB 
b A = 



d<p A 
df 

d<p A 



H — bip J 



H-b<p A ) df. 



dip 



B 



H - b<p B ) df 



df 



We denote the linear equations in the matrix form 
Kx = b 



(70) 



with K 



( fell fel2 
fe21 fe22 



feln \ 
fe2n 



/ n 1 \ 
n 2 



\ fenl fen2 • • • fenn / 



and B 



( h \ 
\b N J 



The system of linear equations is solved by the Singular Value Decomposition (SVD) method [T9] 
due to the singularity of the stiffness matrix. The stiffness matrix is decomposed as a product 
of three matrices, i.e. 



K = U ■ W ■ V T = U ■ [diag( Wj )} ■ V 1 



(71) 



where U and V are orthogonal, and W is diagonal. For the singular matrix K, one or more w'jS 
will be zero. If more w'^s are zero, the matrix is more singular. The matrices computed from 
SVD simply give the orthonormal bases of the range and nullspace of the matrix K. Multiple 
solutions will be found because of the singularity of the stiffness matrix. Any solution to ( 69 ) 
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defines an approximate solution to (67). We optimize over these solutions to satisfy the condi- 
tions mentioned below. 



4.2 PROBLEM DESCRIPTION 

Two dynamical systems are examined in our preliminary tests. One is the Lorenz system |17j ; 
the other is the Hald problem 

We recall the ODEs of the Lorenz system: 



x = a(y — x) 
y = 72 — y — xz 
z = xy — (3z 

a = 10, = 8/3, 7 = 25. 
The governing equation to solve for the coarse function takes the form 

an , x du . . on . 

— ■ a(y - x) + — ■ (jx - y - xz) + — ■ (xy - /3z) = a + ML 



(72) 



(73) 



Figure 14 shows a level set of fine variables corresponding to an identical coarse value c. In 
this case, c equals -0.0081. Also, we can see from the figure that there are multiple fine states 
corresponding to one single coarse state. 

Lorenze system 




c = -0.0081 



y x 

Figure 14: A level set of fine variables 



The second problem is a 4D Hamiltonian system proposed by Chorin et al The ODEs of 
this system are 



xi = x 2 

x 2 = -xi (1 + x%) 

24 = -2 3 (1 + £?), 



(74) 



where (21, 22) and (23, 24) represent the displacements and velocities of two nonlinearly coupled 
oscillators. Similarly, the PDE of the Hald problem is of the form 



^•2 2 + — .(-2,(1 + 2^)) + — .24 + ^--(-23(l + 2?) ) = a + 6n. 



(75) 
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The ultimate goal of choosing coarse variable is to generate autonomous coarse response. Notice 
that the mathematical statement we present is a necessary but not sufficient condition towards 



this goal. Moreover, one can hardly find an analytical solution to (67), and the numerical 



approximations do not have C 1 continuity. So the coarse variable needs to be carefully selected 



from the number of solutions to ( 69 ) so that we obtain desirable characteristics in coarse response 
'close' to being autonomous. This selection is based on three criteria: 

Criterion (1). Consistency of coarse trajectories associated with symmetries of the 
fine system: This criterion means the following: suppose the fine system is such that there 
exists a non-empty set C of trajectories which is endowed with a non-empty (for simplicity, 
countable) set of mappings 

Y C = {s C ■ scifHt)) = s c (f(t)) V t iff 1 and f belong to C}. (76) 

It is possible that there could be more than one such set C for a given fine system. We require 
that our computational procedure is successful in 'blindly' capturing as many such functions sc 
as possible, i.e. without explicit knowledge of 'symmetry' class C built into the calculation of 
sc- 

Criterion (2). Satisfactory match of coarse evolutions from distinct fine initial condi- 
tions with constraints: Since the coarse variable is defined as a function of the fine variables, 
there exist multiple fine states associated with one coarse state and its time derivative. We 
deem a II as acceptable if for any two fine trajectories f 1 and f 2 that satisfy 

no /!(()) = no / 2 (o) 

|(no/i)(o) = |(nof)(o). (77) 

The difference of II o J 1 and Ho f 2 is much smaller, in a suitable norm defined later, than the 
difference of J 1 and f 2 . 

Criterion (3). Local convergence of coarse evolutions over a short time: It is natural 



to require that n be such that whenever (77) is satisfied, Ho f 1 and Ho f 2 agree well over short 
periods of time. 



As discussed in the previous section, (70) has numerous solutions since any vector in the 
nullspace can be added to a solution x. Therefore, the solution is a linear combination of two 
parts: the particular part and the homogeneous part, i.e. 

X — Xp ar t -\- OiXfi 0mi (^) 



where ct is an arbitrary scalar. Xp ar t denotes any particular solution to fl70h, and Xh om is any 
vector in the nullspace. The smoothness of the coarse function can be estimated by comparing 
the ratio of the variance to the mean of x par t and Xhom, where we find that x par t generally 
fluctuates more than x^ om . Thus, by changing the value of a we can 'adjust' the smoothness 
of the coarse function. If a is larger, the coarse function is smoother. If a goes to infinity, 
the coarse function is close to a constant. We will not choose an infinitely large a since a 
constant coarse function is trivial. However, by modifying a we can obtain a coarse function 
with desirable smoothness so as to fulfill the requirement of Criteria (2) and (3). 

A measure is introduced to estimate the difference between a pair of trajectories starting 
from distinct fine initial conditions, as illustrated in the equations below 



II f 1 - f 2 ll 

Af := IIJ J 11 - (79a) 

i/s-dl/^l + ll/ 2 !!) 



I c 1 — c 2 1 



Ac := " ". -, (79b) 

1/2 •(||c l || + ||c»||) 
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where 



N —If 

Eco? and (-) : =t/ 



{.)(t)dt. 



(80) 



T denotes the total simulation time. Here f 1 and f 2 refer to the fine trajectories starting from 
two distinct fine initial conditions, respectively. Similarly, c 1 and c 2 present the corresponding 
coarse trajectories, i.e. c l {t) = H(f l (t)). This measure gives a sense of how much the trajectories 
differ from each other. The next section makes use of these measures to discuss the simulated 
results. 

4.3 NUMERICAL TESTS 

We test the properties of the computed coarse variables and present numerical results in this 
section. 



(1). Consistency of coarse trajectories associated with symmetries of the fine sys- 
tem: One symmetry of the Lorenz system is 



x 2 (0) = 


-x^O) 




x 2 (t) = 


-x 1 ^) 


y 2 (o) = 


-2/(0) 


— > 


y 2 (t) = 


-yHt) 


z 2 (0) = 


^(0) 




z 2 (t) = 


z^t). 



(81) 



Figures [X5|a),[T5|b) and[Xo^c) show the fine evolutions. The trajectories with distinct fine initial 
values are indicated with different colors. The rule applies to all following plots. For x and y 
variables, the trajectories are symmetric about x = and y = 0, respectively, while the evolution 



of z is identical in both cases. Figure 15 ^d) shows the corresponding coarse evolutions, where 



the trajectories totally overlap each other, c is the defined coarse variable n(/) as indicated 



in ( 68 ) . It is remarkable that a discretely computed quantity preserves the symmetry property 
described in Criterion (1). Our procedure may be considered a generalization of the idea of 
finding invariants for autonomous ODE systems as, e.g., those considered in [B]. 

Suppose we choose a trajectory (x (t),y (t),z (t)). Then (81) and any one (or all) of the 
following functions sc(f) '■= \x\, or sc(f) '■= \y\, or sc '■= {\x\, \y\, z} T defines a class C 
and the associated set of mappings Y being discussed here. With these definitions in hand, we 
observe from Figure 15 ^d) that II(/ 1 (t)) = Tl(f 2 (t)) V t for f 1 and f 2 belonging to C. 

For the Hald problem, the symmetry has more than one form 



x 2 (0) = 


4(0) 




x 2 {t) = 


x\(t) 


x 2 (0) = 


4(o) 


— > 


x 2 (t) = 


x\(t) 


4(0) = 


-4(o) 


4(t) = 


-4(t) 


4(o) = 


-4(o) 




x 2 (t) = - 


-4(t), 



(82a) 



zf(0) = 


-4(o) 




xl(t) = 


-x\{t) 


4(o) = 


-4(o) 


— > 


4(t) = 


-xl{t) 


4(o) = 


-4(o) 


4(t) = 


-4(t) 


4(o) = 


-4(o) 




x\{t) = 


-4(t), 



(82b) 



4(o) = 


-4(o) 




xj(t) = 


-x\{t) 


4(o) = 


-4(o) 


— > 


4(t) = 


-x\(t) 


4(o) = 


4(o) 


4(t) = 


x\{t) 


4(o) = 


4(o) 




x\(t) = 


x\(t). 



(82c) 
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Once the fine initial values satisfy any of the possible symmetry forms, the fine evolutions 
will maintain the same kind of symmetry, as shown in Figure [16} In this case, if we choose 



sc(f) ■= \x2\V, or s c (f) ■= {\x 3 \, M} T , or s c (f) ■= \x 2 \, \x 3 \, |x 4 |} T , these 

definitions of sc and classes of trajectories generated from any one of (82a) 
to define autonomous coarse functions that we desire. 



(82b), (82c 



serve 
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Figure 16: Test symmetry on the Hald problem (a) Xi(t); (b) x 2 (t); (c) Xs(t); (d) x 4 (i); (e) c(t) 



(2). Coarse evolutions corresponding to distinct fine initial conditions with con- 
straints: We can see from Figure 17 that for fine trajectories that satisfy the constraints (77), 
coarse evolutions match well, even though fine trajectories show a large divergence. By com- 
paring the fine and coarse differences with the same measure, we observe a significantly larger 
fine difference. The ratio Ac/A/ is equal to 0.1306/0.7449=0.1753, where A/ and Ac are 



defined in (79i) and (|79p), respectively. The same constraints are imposed on the fine ini- 
tial data on the Hald problem. Similar results are found as shown in Figure 18, in which 
Ac/A/ = 0.2501/1.1310 = 0.2211. 
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Figure 18: Test coarse evolution with constraints on the Hald problem (a) Xi(t); (b) X2(t); (c) xs(t); 
(d) x 4 {t); (e) c(t) 



(3). Local convergence of coarse evolutions over a short time: Figure 19 shows four 
examples of the coarse evolutions within a shorter time. The constraints (77) are imposed as 
before. We compute the coarse difference till t = 0.3. The coarse differences Ac are 0.0163, 
0.0154, 0.1030 and 0.0765, respectively, which are generally smaller than the coarse differences 
for longer period of time shown in Figure [i~7|^d). However, for any two arbitrary fine initial con- 
ditions without these constraints, the coarse trajectories do not show local convergence. The 
coarse differences equal 0.2136, 0.2184, 0.2325 and 0.2224, respectively, which are more than 10 
times bigger, as indicated in Figure 20 - this is simply to verify that this result is not simply an 



artifact induced by choosing a coarse variable that is a function on the fine phase space with 
negligible variation about a constant value. Similarly, Figure 21 shows the results of the Hald 
problem. We compute the coarse difference till t = 0.4. For two coarse trajectories with the 
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same initial value and its time derivative, the initial difference is considerably smaller compared 
with that from arbitrary fine initial conditions. 





Figure 20: c(t) with arbitrary fine initial conditions on the Lorenz system (a) Case I: Ac = 0.2136; (b) 
Case II: Ac = 0.2184; (c) Case III: Ac = 0.2325; (d) Case IV: Ac = 0.2224 
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Figure 22: c(t) with arbitrary fine initial conditions on the Hald problem (a) Case I: Ac = 0.3019; (b) 
Case II: Ac = 0.2097; (c) Case III: Ac = 0.2929; (d) Case IV: Ac = 0.1881 
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5. Discussion 



We summarize salient features of our work: 

• Our main contribution is the sound demonstration of running-time averages as appropriate 
coarse variables of fast fine-scale dynamics. No other multi-scale computational method in 
the literature ( "Equation- free" approach [TB], Heterogenous Multi-scale Method (HMM) |25[ 
I12j . Young measure theory |6J) have discussed such a possibility. For example, the averaging 
of the discrete KDV-Burgers type equation in [6] required the explicit knowledge of the 
iV/2+1 invariants of the fast flow with period ./V; such information, unfortunately, is generally 
not available for dynamical systems. 

• We have explored three options for constructing coarse evolution equations. Each has its 
own pros and cons: 

— In problems where time-averaged coarse response is insensitive to fine initial conditions, 
PLIM seems to perform reasonably. The timescale of the coarse dynamics constructed 
using PLIM is the averaging period r; while using the Tikhonov approach and Practical 
Time Averaging, we seek the average of the limit dynamics at the timescale 7j, which is the 
timescale of loading period and is usually provided by the physical problem. Suppose we 
have a MD system with atoms vibrating at timescale t and an external loading applied to 
the system with loading period Tj, with measurements of different resolutions r. Normally 
we have t < r < Tj. PLIM allows us to 'observe' the dynamics at different scales. On 
the other hand, in the study of the singularly perturbed system with e — >• 0, the dynamics 
at the intermediate scale is not well defined. However, computation of the invariance 
PDE is a non-trivial task, especially when e — > 0. The computational methods we have 
tried are quite simplistic and it remains to be seen how Sacker's |21j elliptic regularization 
existence guarantees fares for PLIM with a small number of coarse variables as well as more 
sophisticated numerical schemes tailored to hyperbolic PDE [15] - the latter for the case 
when e is not far less than 1. 

— For the limiting case when e — > 0, the Tikhonov approach is efficient in computation, 
compared with solving the PDE in PLIM and provides a candidate ODE system (defined 
on the slow time scale with no e) for PLIM to be further implemented. However, invertibility 



of the matrix dH/df in (32) is a serious issue. It also has its limitations in that, at best, 
it can approximate time averaged dynamics robustly only when the fast dynamics settles 
on an equilibria. Moreover, even with latitude, it cannot be expected to approximate the 
time-averaged behavior of nonlinear functions of state (while, as demonstrated herein, it 
can on occasion do a reasonable job for linear state functions). 

- Practical Time Averaging takes into account more general cases with the only requirement 
that the fast flow (for each fixed coarse state) remain in a bounded region of phase space. 
We have demonstrated several such cases in this paper. It has been shown that the limit be- 
havior of nonlinear functions can be approximated as well. It is also efficient in computation 
since it allows a very big interval T as the coarse time step. 

Construction of autonomous coarse system is still an open question. With the implemen- 
tation of PLIM, the emergence of physically history-dependent behavior from underlying 
current-state-dependent behavior is observed (for instance, in |22j), which arises as a direct 
consequence of reduction or coarse observation of the fine system. As discussed in detail 
in [23] such dependence is not expected to go away even when time-averaged coarse variables 
are used. For non-history dependent coarse variables, such non-autonomous coarse behavior 
was most explicitly noted in [8] and subsequently in [7]. The computations in [6] do not 
emphasize such dependence but, as pointed out in [23J, it is of course there and is accounted 
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for by the special choice of initial conditions that is made to computationally generate the 
Young measure averages. This issue necessarily exists in HMM too in their reconstruction 
step. For instance, to obtain the stress (the only unknown data) in the macroscopic model 
of complex fluids (in which the coarse variables are the macroscopic velocity field and the 
stress), they |20j apply the Irving-Kirkwood formula, which is dependent on the positions 
of particles. These quantities are defined in the microscopic model. From a physical point 
of view, fine quantities that correspond to a unique macroscopic velocity field can not be 
ensured to generate a unique stress. 

• In the second part of the paper, we consider the more classical question of coarse variables 
that are 'instantaneous' functions of the fine variables. An interesting finding here is that 
we are able to numerically define some such functions whose evolution is guaranteed to be 
unique when evaluated on whole classes of fine trajectories. 
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